arXiv:1509.03655vl [math.AP] 9 Sep 2015 


A COUPLED SURFACE-CAHN-HILLIARD BULK-DIFFUSION SYSTEM 
MODELING LIPID RAFT FORMATION IN CELL MEMBRANES 

HARALD GARCKE, JOHANNES KAMPMANN, ANDREAS RATZ, AND MATTHIAS ROGER 


Abstract. We propose and investigate a model for lipid raft formation and dynamics in biolog¬ 
ical membranes. The model describes the lipid composition of the membrane and an interaction 
with cholesterol. To account for cholesterol exchange between cytosol and cell membrane we cou¬ 
ple a bulk-diffusion to an evolution equation on the membrane. The latter describes a relaxation 
dynamics for an energy taking lipid-phase separation and lipid-cholesterol interaction energy into 
account. It takes the form of an (extended) Gahn-Hilliard equation. Different laws for the ex¬ 
change term represent equilibrium and non-equilibrium models. We present a thermodynamic 
justification, analyze the respective qualitative behavior and derive asymptotic reductions of the 
model. In particular we present a formal asymptotic expansion near the sharp interface limit, 
where the membrane is separated into two pure phases of saturated and unsaturated lipids, re¬ 
spectively. Einally we perform numerical simulations and investigate the long-time behavior of 
the model and its parameter dependence. Both the mathematical analysis and the numerical 
simulations show the emergence of raft-like structures in the non-equilibrium case whereas in the 
equilibrium case only macrodomains survive in the long-time evolution. 


1. Introduction 

Phase separation processes that lead to microdomains of a well-defined length-scale below 
the system size arise in various physical and biological systems. A prominent example is the 
microphase separation in block-copolymers [7] or other soft materials, characterized by a fluid- 
like disorder on molecular scales and high degree of order at larger scales. Here micro-scale 
pattern arise by a competition between thermodynamic forces that drive (macro-)phase separation 
and entropic forces that limit phase separation. Different mathematical models for microphase 
separation in di-block copolymers have been developed, in particular built on self-consistent 
mean held theory (see for example |35| and the references therein) or density functional theory 
developed in [SaiMlE] that leads to the so-called Ohta-Kawasaki free energy. Diblock-copolymer 
type models have also been studied on spheres and more general curved surfaces [321 ED [ini EH 
and show a variety of different stripe- or spot-like patterns. In contrast to materials science 
microphase separation on biological membranes is much less understood, in particular in living 
cells. In this contribution we present and analyze a model for so-called lipid rafts, that represent 
microdomains of specific lipid compositions. 

The outer (plasma) membrane of a biological cell consists of a bilayer formed by several sorts 
of lipid molecules and contains various other molecules like proteins and cholesterol. Besides 
being the physical boundary of the cell the plasma membrane also plays an active part in the 
functioning of the cell. Many studies over the past decades have shown that the structure of 
the outer membrane is heterogeneous with several microdomains of different lipid and/or protein 
composition. Such domains have been clearly observed in artificial membranes such as giant 
unilamellar vesicles gHiiin]. Here a dess fluidic’ (liquid-ordered) phase of saturated lipids and 
cholesterol separates from a ‘more fluidic’ (liquid-disordered) phase of unsaturated lipids. The 
nucleation of dispersed microdomains is followed by a classical coarsening process that leads to 
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coexisting domains with a length scale of the order of the system size. The situation is much 
more complex and much less understood for living cells. Their plasma membrane represents 
a heterogeneous structure with a complex and dynamic lipidic organization. The formation, 
maintenance and dynamics of intermediate-sized domains (10 - 200 nm on cells of jim size j30| ) 
called dipid rafts’ are of prime interest. These are characterized as liquid-ordered phases that 
consist of saturated lipids and are enriched of cholesterol and various proteins miiH]. These rafts 
contribute to various biological processes including signal transduction, membrane trafficking, 
and protein sorting m- It is therefore an interesting question to study the underlying process, 
by which these rafts are generated and maintained, and to understand the mechanism that allows 
for a dynamic distribution of intermediate-sized domains. 

Several phenomenological mesoscale models have been proposed, see for example the review 
m- One class of models argues that raft formation is a result of a thermodynamic equilibrium 
process. Here one contribution is a phase separation energy that would induce a reduction of 
interfacial size between the raft and non-raft phases. The observation of nano-scale structures 
is then explained by including different additional energy contributions. One proposal is that 
thermal fluctuations near the critical temperature for the phase separation are responsible for 
the intermediate-sized structures. Other explanations consider interactions between lipids and 
membrane proteins that could act as a kind of surfactant or could ‘pin’ the interfaces due to their 
immobility [SH uni ES] • Finally, raft formation could also be stabilized by induced changes of the 
membrane geometry and elasticity effects. However, as argued in CHI, all such models are not able 
to reproduce key characteristics of raft dynamics; non-equilibrium effects essentially contribute to 
raft formation. Of particular importance are active transport processes of raft components that 
allow to maintain a non-equilibrium composition of the membrane. The competition between 
phase separation and recycling of raft components is argued to be of major importance for the 
dynamics and structure of lipid rafts. 

Foret m proposed a simple mechanism of raft formation in a two-component fluid membrane. 
This model includes a constant exchange of lipids between the membrane and a lipid reservoir 
as well as a typical phase separation energy. While relaxation of the latter tends to create large 
domains, the constant insertion and extraction of lipids in the membrane ensures indeed the 
formation of rafts. The emerging microdomains are static (in contrast to the lipid rafts on actual 
cell membranes) and the size distribution of these rafts is rather uniform, whereas in vivo cell 
membranes show a dynamic distribution of rafts of different sizes (see the concluding remarks in 
1201 and the discussion of Forets model in USD- Moreover, as already noticed in [18], whereas 
Forets model is motivated by including non-equilibrium effects his model can be equivalently 
characterized as relaxation dynamics for an effective energy given by the Ohta-Kawasaki energy 
of block-copolymers |39| that is well-known to generate phase-separation in intermediate-sized 
structures. 

A similar model by Gomez, Sagues and Reigada [22| considers a ternary mixture of saturated 
and unsaturated lipids together with cholesterol and studies the interplay of lipid phase sepa¬ 
ration and a continuous recycling of cholesterol. An energy that is determined by the relative 
concentration (j) of saturated lipids and the relative cholesterol concentration c is proposed that 
in particular includes a phase separation energy of Ginzburg-Landau type for the lipid phases 
and a preference for cholesterol-saturated lipid interactions over cholesterol-unsaturated lipid 
interactions. The dynamics include an exchange term for cholesterol that is given by an in-/out- 
flux proportional to the difference from a constant equilibrium concentration of cholesterol at the 
membrane. 

Our aim here is to propose an extended model and to present both a mathematical analysis and 
numerical simulations. Similar as in m one ingredient of our model are energetic contributions 
from lipid phase separation and lipid-cholesterol interaction. In addition, we include the dynam¬ 
ics of cholesterol inside the cytosol (the liquid matter inside the cell) and prescribe a detailed 
coupling between the processes in the cell and on the cell membrane. In particular, the outflow 
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of cholesterol from the cytosol appears as a source-term in the membrane-cholesterol dynamic 
and will be characterized by a constitutive relation. We will investigate different choices of this 
relation and will illustrate the implications on the emergence of microdomains. 


1.1. A lipid raft model including cholesterol exchange and cytosolic diffusion. To give 
a detailed description of our model let us fix an open bounded set B with smooth boundary 
r = dB representing the cell volume and cell membrane, respectively. Let (p denote a rescaled 
relative concentration of saturated lipid molecules on the membrane, with cp — 1 and (p = —1 
representing the pure saturated-lipid and pure unsaturated-lipid phases, respectively. Moreover 
let V denote the relative concentration of membrane-bound cholesterol, where v = 1 indicates 
maximal saturation, and let u denote the relative concentration of cytosolic cholesterol. We then 
prescribe a phase-separation and interaction energy of the form 

(1.1) ^ (||Vr¥^P + + ^i2v - 1 - 


with W a double-well potential that we choose as W{ip) — ^ constants 6:, 5 > 0, and 

Vr denoting the surface gradient. The first two terms represent a classical Ginzburg-Landau 
phase separation energy, whereas the third term models a preferential binding of cholesterol to 
the lipid-saturated phase. We define the chemical potentials 


(1.2) 

5F 

^ 5ip 

(1.3) 



—e/S.Y<p + e ^W\(p) — 5 ^{2v — l — (p)^ 

^(2t - 1 -(/p). 


where Ap denotes the Laplace-Beltrami operator on T, and prescribe for the dynamics of the 
concentrations over a time interval of observation (0, T) the following system of equations 


(1.4) 

dtu = D/!\u 

in B X (0,T], 

(1.5) 

—DVu ■ u = q 

on r X (0,T], 

(1.6) 

dtp = ArM 

on r X (0,T], 

(1.7) 

4 2 

dfV = ArO + q = ^Apv - ^Arp + q 

0 0 

on r X (0,T], 


where u denotes the outer unit-normal field of S on T. The system is complemented with 
initial conditions foi* R? ^ ^^5 respectively. The first equation represents a simple 

diffusion equation for the cholesterol in the bulk, the second equation characterizes the outflow 
of cholesterol. The third equation is a Cahn-Hilliard dynamics for the lipid concentration on the 
membrane, whereas the equation for the cholesterol on the membrane combines a mass-preserving 
relaxation of the interaction energy and an exchange with the bulk reservoir of cholesterol given 
by the flux from the cytosol. This combination yields a diffusion equation with cross-diffusion 
contributions and a source term. 

To close the system it remains to characterize the exchange term q. We follow here two 
possibilities: First we prescribe a constitutive relation by considering the membrane attachment as 
an elementary Reaction’ between free sites on the membrane and cholesterol, and the detachment 
as proportional to the membrane cholesterol concentration, expressed by the choice 


(1.8) 


q — ciu{l — v) — C 2 V. 


A similar coupling of bulk-surface equations has been investigated in a reaction-diffusion model 
for signaling networks [mils]. As a second possibility we consider choices of q that allow for a 
global free energy inequality for the coupled membrane/cytosol system. These two different cases 
could be considered as a distinction between open and closed systems and one important aspect 
of this work is to evaluate the consequences of these choices for the formation of complex phases. 
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We remark that the system conserves both the total cholesterol and the lipid concentrations 
since for arbitrary choices of q we obtain the relations 


(1.9) 

( 1 . 10 ) 


UL 


udx + J V dl-L^ = 0, 


Let us contrast the above model with the Ohta-Kawasaki model for phase separation in diblock 
copolymers mentioned above. Let Q C denote a spatial domain, (p the relative concentration 
of one of the two polymers and let m j- (p he the prescribed average of Lp over ft. The mean 
field potential z is then given by 


-Az = p) 


m 


in f], 


Vz • Po = 0 on 


I 

Jq 


z = 0. 


Then a free energy is prescribed of the form 

( 1 . 11 ) = 






dx + 




z\‘^ dx, 


where a > 0 is a fixed constant. Note that the last term can also be written as ||| 
A relaxation dynamics of Cahn-Hilliard type is then typically considered given as 

6Tok 




( 1 . 12 ) 


dtif = An, 


= 


6(p 


= —£Ap>-\-e crz 


The energies T and Tqk both contain a Cahn-Hilliard type energy contribution that favors 
macro-phase separation but are different in the additional terms. However we will see below that 
stationary patterns for our lipid raft model with the choice (1.8) are for small d > 0 closely related 
to stationary points of Tqk- If on the other hand one considers simple choices for q that lead to 
a global free energy inequality, we will observe a macro-scale separation of all saturated lipids in 
one connected domain. This indicates that in fact non-equilibrium processes are responsible for 
lipid raft formation. 

Several asymptotic regimes are interesting in view of the different parameters included in our 
model. We will investigate the limit e ^ 0 that corresponds to a strong segregation limit and 
leads to a model where no mixing of the lipid phases is allowed and the domain splits into regions 
where p) = 1 and p? = —1 respectively. This limit corresponds to the sharp interface limit in 
phase-field models and connects to the analysis of the Ohta-Kawasaki model in [3Hj. Since the 
cytosolic diffusion in biological cells is known to be much faster than lateral diffusion on the cell 
membrane another natural reduction of the model appears in the limit D ^ oo that leads to a 
non-local model defined solely on the cell membrane. Finally, assuming the effect of the lipid 
interaction with cholesterol to be large motivates to consider the asymptotic regime 6 I 0. 


1.2. Outline of the paper and main results. In the next section we will derive the model 
from thermodynamic considerations. In particular we will show that for arbitrary 
choices of the exchange term q the surface equations and the bulk (cytosolic) equations are 
thermodynamically consistent when viewed as separate systems. Depending on the specific choice 
of q we may or may not have a global (that is with respect to the full model) free energy inequality. 
We will present examples for both cases. 

In Section we will first derive a reduced raft model in the large cytosolic diffusion limit by 
formally taking D ^ oo. We then analyze the qualitative behavior of the (reduced) system in 
terms of a characterization of stationary points and an investigation of their relation to stationary 
points of the Ohta-Kawasaki model. A formal asymptotic expansion for the sharp interface 
reduction s ^ 0 of our raft model is presented in Section Here we also briefly discuss the 
resulting limit problem that takes the form of a free boundary problem of Mullins-Sekerka type 
on the membrane with an additional coupling to a diffusion process in the bulk and including 
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an interaction with the cholesterol concentration. In Section we present numerical simulations 
of the full and reduced raft model. In particular, we study spinodal decomposition, coarsening 
scenarios and the possible appearance of raft-like structures as (almost) stationary states for 
different choices of the exchange term q and for different parameter regimes. Some conclusions 
are stated in the final section. 

2. Thermodynamic justification of the lipid raft model 

In this section we will derive the governing equations for the lipid raft model from basic 
thermodynamical conservation laws using a free energy inequality. Our arguments are similar to 
an approach used by Gurtin 121ES], who derived the Cahn-Hilliard equation in the context of 
non-equilibrium thermodynamics. We first of all consider the equations which have to hold on 
the surface, will then consider the equations in the bulk and subsequently we will couple both 
systems. 

The basic quantities on the membrane surface are the rescaled relative concentration of the 
saturated lipid molecules (^, the concentration v of the membrane-bound cholesterol, the mass 
flux of the lipid molecules, the mass flux Jy of the cholesterol, the mass supply of cholesterol g, 
the surface free energy density / and the chemical potential fi related to the lipid molecules and 
the chemical potential 9 related to the surface cholesterol. The underlying laws for any surface 
subdomain S C F are the mass balance for the lipids 

( 2 . 1 ) ^ [ (pd'H‘^ = - f J^-ndT-L^, 

the mass balance for the surface cholesterol 

( 2 . 2 ) ^ f vdH^^— f Jy-ndV}^ f qdH^ 

dt Jy: Jas JS 

and the second law of thermodynamics, which in the isothermal situation has the form 

(2.3) 4 f fdV^ < - [ ■ n) + OJv ■ n) dV} + f Oqdn^ 

dt Jy Jot. js 

where / denotes the surface free energy density and where n is the outer unit conormal to dY^ 
in the tangent space of F. In addition, we denote by dJ-L^ the integration with respect to the 
d—dimensional surface measure and /,Vrv^ denotes the partial derivatives of / with respect to the 
variables related to in a constitutive relation / = /(..., VrV^, .••)• Similarly we will denote 
with a subscript comma partial derivatives with respect to other variables. For a discussion of 
these laws in cases in which source terms are present and at the same time / does not depend on 
Vr^ we refer to [24] . [26l Chapter 62], and [40]. Thermodynamical models of phase transitions 
with an order parameter typically involve a free energy density / which does depend on Vr^- 
In this case the free energy flux does not only involve the classical terms and 9Jy but also 
a term This is discussed in [25l O [2]. Gurtin [25] introduces a microforce balance 

involving a microstress ^ in order to derive the Cahn-Hilliard equation. Here, we do not discuss 
the microforce balance and instead already use the form ^ which could be derived in 

our context in the same way as in [25]. However, in order to shorten the presentation we do not 
state the details. We hence obtain ( |2.3[ ) as the relevant free energy inequality in cases where no 
external microforces are present. 

Since the above (in-)equalities (I 2 IJ-Q hold for all S, we obtain with the help of the GauB 
theorem on surfaces the local forms, compare [25], 


(2.4) 

dfif + divr = 0 

in F X (0,T], 

(2.5) 

dtv + divr = q 

in F X (0,T], 

(2.6) 

dtf + divr(M-C - + OJv) < Oq 

in F X (0,r]. 


6 


H. GARCKE, J. KAMPMANN, A. RATZ, AND M. ROGER 


With the constitutive relation / = /(p, (/p, VrV^) we obtain from the local form of the free energy 
inequality 


f,cpdt^ + /,'i;5tP+Vr/i •«/(/? + Vr^ • Jv~\~ 

(divr + (divr - dt (^divr/,Vr^ ^ 
Using the conservation laws ( |2.4| ) and ( |2.5[ ) we obtain 

if,if - divr(/,Vrv^) ~ + {f^v - 0)dtv + Vr/i • J^p + Vr^ ' Jv 


The fact that solutions of the conservation laws with arbitrary values for dt^p and dtv can appear 
is used in the theory of rational thermodynamics to show that the factors multiplying dt^p and dtv 
have to disappear as they do not depend on dtp> and dtv. We refer to Liu’s method of Lagrange 
multipliers [3lj and to for a more precise discussion on how the free energy inequality 

can be used to restrict possible constitutive relations. We now choose the following constitutive 
relations which guarantee that the free energy inequality is fulfilled for arbitrary solutions of (2.4), 
In fact, choosing 


(2.7) 

(2.8) 

(2.9) 

( 2 . 10 ) 


/^ = A - divr(/,vrv^)> 

0 = f,v, 

Jf) — 

Jy = —DyS/^O 


with Dy >0 ensures that the free energy inequality is fulfilled for all solutions of the conser¬ 
vation laws. More general models, e.g. taking cross diffusion into account, are possible and we 
refer to [5] for an approach which can be used to obtain more general models. 

We now consider the governing physical laws in the bulk. As variables we choose u which is 
the relative concentration of the cytosolic cholesterol, the bulk chemical potential the bulk 
free energy density fb{u)^ the bulk flux Jy and the surface mass source term Qy. We need to fulfill 
the following mass balance equation in integral form which has to hold for all open U C B: 

( 2 . 11 ) ^ f udx^— f Jy'udT-L^-l- [ qudV? 

ht Ju J{du)\r J{du)nr 

and the free energy inequality 

(2.12) ^ f fb{u)dx<- f /i„J„ • + f iiuqudT-L^ 

Ju J(du)\r J{du)r\r 


which also has to hold for all open U C B. Here we allowed for a source term qy on L and in 
accordance to our discussion above we introduced the free energy source jiyqy in (2.12). As on 
the interface the free energy source term is given classically as a product of the mass source term 
and the chemical potential. If we use the fact that we can choose an arbitrary open set U which 
is compactly supported in B (which then implies {dU) H L = 0) we obtain with the help of the 
GauB theorem 


(2.13) dtu-\-diYjy — Q inSx(0,T], 

(2.14) dtfb{u) div {fly Jy) <0 in S x (0,r]. 

Choosing 

(2.15) i^u = fUu), 

(2.16) Ju = -M{u)V{fl{u)) 

makes sure that (2.14) is true for all solutions of ( |2.13D . In the following we will often choose 
M{u) — and this will lead to the linear diffusion equation ( |1.4| ) 

dfU — DyAu = 0 in H X (0,T]. 








SURFACE CAHN-HILLIARD MODELING RAFT FORMATION 


7 


Choosing U such that dU H F ^ 0 we obtain from ( 2.11[ ) 

0= / {dfU-\-div Ju) dx = / {Qu ^ Ju • i^)dT-L^ 

Ju Jdunr 

which gives, since U is arbitrary, 

Qu = - Ju • u (= DyVu • u) in r X (0,r], 

and which yields (1.5) for q = —Qu- We will now state global balance laws in the case where the 
mass supply for the interface stems from the bulk and vice versa. 

Lemma 2.1. We assume that the above stated mass balanee equations hold for the bulk and the 
surfaee and assume in addition that qu — —q- Then it holds 




= 0 — ^dH^ = 0. 

at Jr 


« jfi 

If in addition, the free energy inequalities in the bulk and on the surfaee are true it also holds 

^{Jr Jb ~ iau)dJ-r. 

Proof The first two equations follow from (2.1), ( |2.2[ ) with S = F and ( |2.11| ) with U = B since 
dU = F and since 5F = 0. The total free energy inequality follows similarly from (2.3) and (2.12). 

□ 


Remark 2.2. (i) In the above lemma we chose qu = —that is the mass lost on the surface 
generates a source of mass for the bulk. 

(ii) Several constitutive laws for q make sense. It is possible to consider 

(2.17) q^-c{9-yiu), c>0 


which leads to model for which the total free energy decreases. In this case we obtain 

Ib L ~ ~ 

(hi) We also consider the reaction type source term, compare ( |1.8[ ), 

q — ciu{l — t) + C 2 V. 


Also in this case we obtain a consistent model which fulfills the bulk and surface free energy 
inequalities with source terms as stated above. However, in this case the total free energy as 
the sum of the bulk and the surface free energy might increase which can be due to the fact 
that we neglect energy contributions generated by the detachment and attachment process, 
(iv) One possible choice for the surface free energy density is, compare (0, 

2 e Id 

For fh{u) = qu = —q and Du — D, ^ — Dy — 1 we obtain the system ( |1.2| )-( [T77l ). 

With the above quadratic choice for and arguing as in the derivation of the energy inequality 
one can prove the identity 


(2.18) 


+ f (iVr^lM [ D|Vn|2= f q{e - u) 

Jr J B Jr 


where the last term is non-positive for the choice (2.17). Any other evolution that is based 
on a choice of q such that the right-hand side of (2.18) is always non-positive decreases the 















H. GARCKE, J. KAMPMANN, A. RATZ, AND M. ROGER 


total free energy. This is in particular the case for choices of q such that doii-P^ can be 
characterized as a gradient flow. 

In the following we are mainly interested in the dependence on the parameters 6:, 5, and 
therefore as in |(iv)| above we always set Du — j = 1, in which case the above 

choices of free energy densities and mass fluxes yield the system (1.2)-(1.7). 

3. Qualitative behavior 

In this section we will investigate qualitative properties of the model and of the 

asymptotic reduction in the large cytosolic diffusion limit that we derive below. One key question 
here is whether or not our lipid raft model supports the formation of mesoscale patterns. We 
will distinguish different choices for the exchange term q and compare evolutions that reduce the 
total free energy with the evolution for the choice q given by the reaction-type law (1.8), which 
we consider as a prototype of a non-equilibrium model. We remark that most of the arguments in 
this section are purely formal; a rigorous justification is out of the scope of the present paper. In 
particular, we assume the existence of smooth solutions, their convergence to stationary states as 
times tends to infinity, and that the long-time behavior of the full system asymptotically agrees 
with that of the reduced system developed below. 

We first observe that under an additional growth assumption on the exchange term q we obtain 
energy bounds, even in the case that the system does not satisfy a global energy inequality. 

Proposition 3.1. Assume that q has at most linear growth, that is there exists A > 0 sueh that 

(3.1) |g((/p,R,T)| < A(1-h |v:^| + |r| + |t|) /or a// (/p,r,t G M. 

Then for all 0 < t < T and all D > > 0, 0 < 6: < Sq any solution of (irii-Q with initial 

data (/pojRQjL’o satisfies 

(3.2) + ^ - C'('^ATTo,£o,t’o,V^o,^^o)- 

Proof From ( |2.18D we deduce that 

(3.3) 

JB Jr 

For the last term we use the estimate 


{9 - u)q{(p,u,v) < J + C/i{l + (p‘^ + + v‘^)j 


(3.4) 


— y* +c'a j r +1) J + c'a j 


For the second term on the right-hand side we obtain 

(3.5) - C{l + J^W{ip)) < C{so){l + T{v{;t),^{;t))). 

Using [28l Chapter 2, (2.25)] the third term on the right-hand side of ( |3.4[ ) can be estimated by 
bulk quantities, 

^ l^\Vuf + C{D) 

for a suitable constant C{D), which in particular yields for all D > Do 

(3.6) < J IVnp + C{Do) u^. 
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To estimate the last integral on the right-hand side of (3.4) we first use Young’s inequality and 


(1.3) to obtain 


+ -9{l + ^(1 + < 29‘^ + ^(1 + 


and further deduce that 


(3.7) 


< ^_l02^cil + W{cp)), 
< C'((5,eo)(l + 


We therefore obtain from (3.3)-(3.7) that 


d 

dt ' 


. t)) + ll^ u{; tf) + I ^ (lVr/ip(-, t) + |Vr 0 p(-, t)) 

< C(5, A, Do, so) u{; tf) + C{6, A). 

By the Gronwall inequality we deduce the claim. 


□ 


Note that for the choice (1.8) of q assumption (3.1) is not satisfied. However, for any modification 
that coincides with that choice on a bounded domain in the r, v plane and that has at most linear 


growth outside the conclusion of Proposition 3.1 holds. For (2.17) or any other choice of q that 
implies a total free energy inequality we obtain an even better estimate, since now the right-hand 


side in (2.18) is non-positive. 


3.1. A reduced model in the limit of large cytosolic diffusion. Since in the application 
to cell biology the bulk (cytosolic) diffusion is much higher than the lateral membrane diffusion 
a reasonable reduction of the model can be expected in the limit D - A oo In the case that the 
exchange term q satisfies assumption (3.1) we deduce by Proposition 3.1 that j^D\Vu\^ is 
bounded uniformly in D. The same conclusion holds in the case of any 


evolution by ( |2.18[ ). Therefore in the formal limit D 
and obtain the (non-local) system of surface PDFs 


:ree-energy decreasing 
oc we conclude that u is spatially constant 


(3.8) 

(3.9) 

(3.10) 


dt^p = ArAt 

/i = —eAr<A> + e~^W'{ip) — 5~^{2v — 1 — ip) 
4 2 

dtv — AyO + q = -r Arf - -r Ar<A> + q{p, u, v) 
d 0 


on r X (0, T], 
on r X (0, T], 

on r X (o,r]. 


This system is complemented by initial conditions for p and v. The cholesterol concentration 
u = u{t) is determined by a mass conservation condition 


(3.11) 


u{t) + / t(-, t) = M, 

Jv 


where M > 0 is the total mass of cholesterol in the system. 

Note that the transformation of the coupled bulk-surface system into a system only defined in 
the surface has the price of introducing a non-local term by the characterization of u through the 


mass constraint. The reduction (3.8)-(3.11) is similar to the reduction to a shadow system for 


2x2 reaction-diffusion systems introduced by Keener see also the discussion in 133. 

In the qualitative analysis below and the numerical simulations we will often restrict ourselves 


to the special choice q of the exchange function given in (1.8). For the reduced model it is then 


possible to compute the evolution of the total mass of u and t, which are related by (3.11). We 
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deduce, cf. ( 2 . 11 ), 


(3.12) ^J^u{t)dx = - J q{(p{-,t),u{t),v{-,t)) d'H'^ = J -ciu{t){l - v){-,t) + C2v{-,t) 


dH^ 


= J{ciu{t) + C2)v{-,t)d'H‘^ - 


'B 


u{t) dx 


= {ciu{t) + C2) ( M - 


u{t)dx\ — / u{t)dx 

iB / 1^1 


'B 




+ Cl 


M- in 


C2 ] / u{t) dx + C2M 


\B\ / JB 

and therefore we see that u{t) remains nonnegative if it was initially nonnegative (which is the 
relevant case) and converges for t ^ 00 to Rqo? which is the positive zero of 


p{z) = -Ciz'^ + (ci 


M - |r| . C2M 


thus 


(3.13) 


Ury 



- -) + 

cP ^Cl|i?| 


Since p(0) > 0 and p(y^) < 0 we also obtain that Jp v{t) remains in [0, M] for all times. 


We remark that if q is given as in (1.8), then assumption (3.1) does not hold. Nevertheless we 


can obtain the reduced model for this particular choice of q if we start with a modified version: 
Replace first q by 

q{u^ v) — ciu — ciri{u)v — C2P, 

where 77 : M ^ M is any smooth, monotone increasing and uniformly bounded function with 
77 (r) — r for \r\ < Now the above arguments apply and we obtain the non-local model 

with exchange term q. By analogous computations as above we then deduce 

d f . . _ ci|r| 

_/„(*)*= -T^ 


j u(t)dx+j u{t)dx)+C2^{M - J u{t)dx), 


which yields 0 < f^u{t)dx < M for all t > 0 if this property holds for the initial data. But 
this implies q{u{t)^v{'^t)) — q{u{t)^v{-^t)) for all t > 0. This justifies to consider in the following 
analysis and in the numerical simulations in Section]^ the exchange term q from (1.8) also for 
the non-local reduction. 


3.2. Stationary points. We are in particular interested in the long-time behavior of solutions 
and will therefore next investigate stationary points of our lipid raft system. For the full system 
( |1.2[ )-( [T?7| ) stationary points Roo, Poo) are characterized by the equations 

0 = Ar/ioo on F, 

(3.14) 0 = Ar 6 >oo + <?(v^oo,^oo,^oo) on F, 

Auoo =0 in S, -DVuoo • p = ^(^ 00 ,'^ 00 ,'^ 00 ) on F. 

This implies that is constant and 

11 9 

( 3 . 15 ) —<£:Ar(^Qo H W (^Poo) — t( 277 oo ~ 1 ~ ^ 00 ) T Moo — 1 “ Moo^ 

so 2 

2 

(3.16) ^(v^ 005 '^00 5'^ 00 ) — ~^A(27;oo ~ 1 ~ V^oo)? 

(3.17) Auoo =0 in B, -DVuoo * p = ^(^ 00 ,'^ 00 ,'^^oo) on F. 
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Alternatively, this system can be characterized by 


(3.18) 


1 


6:Af^OO ^ (v^oo) — Qi^Poo) Mo05 


where Q(^oo ) = \ 0oo is a nonlocal function of Lp^o as Rqc'^oo and hence O^o are determined by 
(poo through ( 3.16 ), ( 3.17 ). 

3.3. The case of an energy-decreasing evolution. Let us in the following first consider the 
case that q is chosen such that the right-hand side in (2.18) is nonpositive, hence the total free 
energy is decreasing. For any stationary point Roo?'^oo? the energy inequality (2.18) yields 
that Poo^Ooo and u^o are constant, in particular by ( |3.15 ) 

1 9 

—eAYPoo~\ - W'{poo) — -h/ioo ^ 1^5 

^(V^OO 5 ROO 5 '^Oo) ~ 0- 

In addition, we fix the total lipid and cholesterol masses as they are for any evolution determined 
by the initial conditions. We therefore prescribe, for M, M\ given. 


M = 


I RqO I '^OO — |^|RoO I 

Jb Jr Jr 

L' 


Afl — I Poo 

In particular, the stationary state p^o coincides with a critical point of the Cahn-Hilliard energy 
subject to a volume constraint. The condition q = 0 provides an additional relation between p^^ 
Uoo and L’oo- Ir fhe case of the exchange law ( 2.17| ) this determines fpVoo and 0oo- 

We can elaborate the connection with stationary points of the Cahn-Hilliard equation a bit 
more if we in addition assume that (uoo^Voo^Poo) is a local minimizer of the energy ( |1.1[ ). We 
represent the latter as 

(3.19) p) = J^i{p) + V^), 


where 


Fi{ip) := + Mv,^) := I^Lp 2 v-l- iff 


We also assume that both Jp p^ and Jp v^o are fixed by the initial data and the condition = 0, 
which holds in particular in case of the exchange law (2.17). 

For any t, p with fpV = Jp v^o and fpP = fp Poo we then have 


J^{2v - 1 - iff = J^(2v - 1 - if - j{2v - 1 - ¥3)) + J (J(2^; - 1 - 

> J (j-{‘^Voc - 1 - <^oo) 

and deduce that under the respective mass constraints the minimizer of given by 

those (v^p) for which 2 t — 1 — (/:? is constant, and that the minimum only depends on and 
Jp p. Since 9oo and thus 2v^ — 1 — p^ are constant, we see that Uoo, Too, ^oo minimizes J^ 2 - 
Moreover, for any p satisfying the mass constraint the minimum of J^2 is attained by = ^(1 + 
p + f( 2 v — 1 — p)) and is independent of p. Therefore poo RS above needs to be a local minimizer 
of the Cahn-Hilliard energy subject to a given mass constraint. In the case of the Cahn- 
Hilliard dynamics in an open convex set in M^, is is known m that stable stationary points 
converge in the sharp interface limit s ^ 0 to configurations with one connected phase boundary 
of constant curvature. In analogy, one therefore might expect that for s > 0 sufficiently small 
the only local minimizer in our lipid raft model for choices q that lead to an energy-decreasing 
evolution are given by configurations with one lipid phase concentrated in a single geodesic ball 
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on r. In particular, such local minimizer do not represent mesoscale pattern like lipid rafts. In 
Section |5.3.4 we present a numerical simulation for energy decreasing dynamics that confirms the 
expected behavior. 


3.4. The case of the exchange term ( |1.8| ). Let us next discuss the choice of q as given in 
(1.8). The representation (3.18) shows some similarity to the equation for stationary points of 
the Ohta-Kawasaki model described above. We can make this more transparent in the case of 
the reduced system (3.8)-(3.11 ), at least if we presume that the long term behavior of the reduced 
systems captures the respective behavior of the full system (which means that the order of limits 
D ^ oo and t ^ oo can be interchanged). Since Rqo is constant we obtain, writing ci = ciu^o for 
convenience, that 


(3.20) 

and further 

^('^005 Pqo) ~ 


0 = 


J q{uoo,Voo) = |r|ci - (ci + C 2 ) J 


Cl + C2 , s Cl + C2 , \ I ~ 

o \^oo 1 ^ 00 ) J-) ¥^ 00 ) “1“ Cl 


S(ci+C2) C 1 +C 2 , 1 ~ 

-- ^00 -—(1 + ‘foo) + Cl 


S{ci + C2) 


5(ci + C 2 ) 


<9oo - 
ffna — 


f 

f 


Ooo 

One 1 — 


<5(ci + C 2 ) 2 


4 

Cl + C2 


(5 

^00 


/<' 


1 ^ 00 ) 


Cl + C2 


(1 + (/Poo) + Cl 


/ 




Cl +C2 


where we have used (3.20) in the third line. This implies by (3.14) 

Since is constant we deduce from equation (3.18) 

—e/S.Y<^oo + c ^W\ipoo) — /ioo 2 2 ~ ^ 

_ \ ^ X n C 1 +C 2 / . 

— Moo 2 / -4— V ~ 

For 6^1 this can be approximated by 
(3.21) 


^00 


/ 


^00 


Cl + C2f ^ , 5(ci+C2)\-1 


V^oo 


/ 


V^oo 


—eAr^Poo + e ^W'{(poo)-\ -;;—(—Ap) ^ 


V^o 


/ 




— Moo + 




The latter equation corresponds to stationary points of the Ohta-Kawasaki functional 


(3.22) 
with 

(3.23) 


•^OA'(V^) = 


2lVv^r + e 




2 

if-i 


a = 


Cl + C2 


where the constant on the right-hand side of (3.21) should be interpreted as a Lagrange multiplier 
associated to a mass constraint for ip. 

The total lipid mass / ^ is given by the initial data. We can also identify ci + C 2 as a function 
of the data. First its value is characterized by Rqo and using (3.13) we deduce 


(3.24) 


Cl +C2 


1 / M |F| 
“ 2V “ \B\ 


, , , Pf M |r| \2 |r| 

-TT^Ci]+. -(c2 + -^Ci--^Ci) +C1C2^. 
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In Section 5.3.2 we will present simulations that confirm that the long-time behavior of the reduced 
system is for 5 <C 1 in fact very close to that of the Ohta-Kawasaki dynamics. In particular, 
in contrast to any choice of q that induces an energy-decreasing evolution, in the case of the 


exchange term (1.8) we in fact see the occurrence of mesoscale patterns. 


Remark 3.2. Let us highlight one key difference in the long-time behavior of our model in the 
different cases considered above. For a free-energy decreasing evolution stationary states are 
characterized by the properties that O^o is constant and q^o — q{^oo^ Rqo? '^oo) zero. For the choice 


(1.8) of the exchange term on the other hand and the reduced system we have stationary states 


with non-vanishing q^^^ which correspond to an persisting exchange of cholesterol between bulk 
and cell membrane. This process eventually allows for the formation of complex pattern on a 
mesoscopic scale. 


4. Sharp interface limit ^ 0 by formally matched asymptotic expansions 


In this section we formally derive the sharp interface limit of the diffuse interface model (1.2)- 
0. We assume throughout this section that the tuple (r^, /i^, solves (1.2)- 

0 to a limit (r, /i, 9). Furthermore, we suppose that the 

0 to a sharp interface. This 


(1.7) as 6: 

(1.7) and converges formally as 

family of the zero level sets of the functions cps converges as 6: 
interface, at times t G [0,T] is supposed to be given as a smooth curve C F which separates 
the regions = 1} and {(p{'^t) = —1}. By a formal asymptotic analysis, we conclude that 

the limit functions (r, (^, t, /i, 9) can be characterized as the solutions of a free boundary problem 
on the surface F, which is again coupled to a diffusion equation in the bulk B. Other examples 
for this method and more details on formal asymptotic analysis can be found in El n El [ig El] 
which is by far not a comprehensive list of references. 

The obtained limit problem describes a time-dependent partition of the surface F into different 
phase regions 

(4.1) F+(t) := = 1} and T~{t) := = -1} 


and the dynamic of the interface 7 (t) between these two regions. We denote the corresponding 
separation of the surface-time domain F x (0,r] by 

U := {(x,t) ev X {0,T] : X e U(t)}. 

The limit problem then takes the following form. 


Sharp Interface Problem 4.1. The sharp interface model obtained from the formal asymptotic 
analysis is given by 


(4.2) 

(P = ±1 

on U, 

(4.3) 

dtu = DAu 

in R X (0,T], 

(4.4) 

—DVu ' v — q ciu{l — v) — C 2 V 

on F X (0,T], 

(4.5) 

o 

II 

< 

on U, 

(4.6) 

dtv = ArO + q 

on U, 

(4.7) 

e=‘^{2v-lTl) 

on r^, 

(4.8) 

2/i 9 = cqKq 

on 7 , 

(4.9) 

o 

II 

+ 1 

on 7 , 

(4.10) 

o 

II 

+ 1 

on 7 , 

(4.11) 

-2V = [VrM]l • 

on 7 , 

(4.12) 

-V = [Vr^jl • lyy 

on 7 , 









14 


H. GARCKE, J. KAMPMANN, A. RATZ, AND M. ROGER 


where is the jump across the interface 7 and z/-^(xo,to) G Tx^V denotes the unit normal to 
7 (^ 0 ) in xo G 7 (^ 0 ), pointing inside r+(to). The geodesic curvature of ^{t) in V is denoted by 
and V(xo,to) denotes the normal velocity of 7 (^ 0 ) in xq G 7 (^ 0 ) in direction of z/-^(xo,to). 
For its precise definition, let • U ^ ^{t) C F, t G (to — 5, to + 5) be a smoothly evolving family 
of local parameterizations of the curves 7 (t) by arc length over an open interval [/ C M and let 
7 to('^o) = ^ 0 - Then the normal velocity in (xo,to) is given by 


see also HB. 


V(x„,(o) = ^ 


7t(so) • 
to 


We first deduce the existence of transition layers between the phase regions. By assuming that 
the functions admit suitable expansions with respect to the parameter s in the 

transitions layers and in the regions away from the interface respectively, we can then deduce 
that the limit functions at least formally need to fulfill (4.2)-(4.12). 


4.1. Asymptotic analysis: Existence of transition layers and outer expansion. We 

start our analysis by expanding the solutions to the coupled model in the outer regions, where ips 
attains values away from zero. We assume that in these regions all functions in ( |1.2[ )-( [T?7| have 
expansions of the form 

00 

fe = Y.^^fk 

where fs = (fe, ^ 6 , • • •, etc. 

Since we postulated the existence of different phase regions characterized by the values of the 
limit function we should first address the existence of these phase regions. To this end, we 
collect all terms of order in 0 and obtain 

fF'((^o) = 0 


and as a consequence we obtain p)o = ±1 as the only stable solutions. Since p^o is the dominant 
term in the assumed expansion as ^ 0 , we deduce ^ ±1 as 6 : ^ 0 and the existence of the 
claimed transition layers. This justifies (4.1) and ( |4.2D . 

The discussion of equations ( |1.3| ) and (1.4)-( [T?7| ) is then straightforward. Comparing the terms 
of order 0(1) in their corresponding equations allows us to deduce (4.3)-(4.7). 

Due to the (possibly steep) transition between the regions F+ and F“ near the interface, the 
spatial derivatives occurring in the system might contribute terms which are not necessarily of 
order 0(1) (with respect to s) in a neighborhood of the interface. This motivates the need for a 
more detailed analysis of the functions in the neighborhood of the interface which we will address 
in Section IT^ 


4.2. Coordinates for a neighborhood of the interface. As stated above, we suppose that the 
zero level sets of converge to some (smooth) curve 7 (t) with inner (wrt. F+(t)) unit normal field 
We then introduce on a small tubular neighborhood N of 7 (t) a new coordinate system 
which is more suitable for the analysis in the transition layer. We remark that the construction 
of these coordinates presented here is more complicated which is due to the fact that A is a 
neighborhood of 7 (t) in the manifold F. For a similar example, we refer the reader to fT7] . 

As above, let 7 ^ : t/ ^ 7 (t) be a local parametrization of the curve ^{t) by arc length over an 
open interval ?7 C M. It is then possible to extend 7 ^ to a local parametrization of N by means 
of the exponential map from differential geometry. 

While details on this map can be found in the literature nacs], it is sufficient for our purpose 
to quickly recall its definition. For a given point p G F and a vector a G T^F, there exists a 
unique geodesic curve such that c^(0) = p and c^(0) = a. The exponential map exp^ in p is 
then defined for all a G TpT for which c^(l) exists and is given by exp^(a) = c^(l). Note that for 
2 : G [0,1] one can easily check exp^( 2 ;a) = 
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The distance between a point x G T and the interface 7 (t) is defined as 
d{x^ t) inf{/(c) |c : / ^ r, c connects x and ^{t) } 
where /(c) denotes the length of the curve c. Setting 

= exp^^(^)(ri/^(74(s),t)), 

we obtain that is a parametrization of a neighborhood N{t) of 7 (t). If we choose N{t) small 
enough, the properties of the exponential map imply that r is the signed distance between the 
point Tt( 5 ,r) and 7(t), that is for {x^t) — Tt( 5 ,r) G N{t) we have 


r = d{x^ t) 


ifrGr+(t) 


For the asymptotic analysis, it is necessary to adapt the parametrization to the length scale of 
the transition layers. We therefore use the rescaled parametrization 

(4.13) A{s,z,t) = 

of N{t)^ where z = ^. In particular, z can be written as a function of x and t. 

Let us remark that as 7 (t) is the zero level set of the signed distance function d, the tangent 
space of 7 (t) is the subspace of the tangent space of F which is orthogonal to the surface gradient 

Vrd and thus the normal vector of 'j{t) is given by . Because of 


0 = —d{'jt{s),t) = Vrd(7t(s),t) • dt'jtis) + dtd{jt{s),t) 


we can thus compute 


dat{s) • u = - 


dtd{'yt{s),t) 


Vrd(7t(5),t) 

and therefore the time derivative dfZ on ^{t) fulfills 

dtz = 

Remark 4.2 (Gradient, Divergence and Laplace Operator in the new coordinates). From the 
definition of A in ( |4J3l ) we see 

(4.14) A{s,0,t) = jt{s), 

(4.15) dsA{s,0,t) = ds'ftis) = 7 ^( 5 ) and 

(4.16) dzA{s,0,t) = eu^i'ytis)). 

Furthermore, the curve r i-A Tt(s,r) is geodesic by definition and hence we have 

VdzzA{s,z,t) = 0 

where V is the projection on the tangent space 7A(s,z,t)r on F in A{s,z,t). 

These observations allow us to calculate 

dz {dsA{s, z, t) ■ dzA{s, z, t)) = dzsA{s, z, t) ■ dzA{s, z, t) + dsA{s, z, t) ■ dzzA{s, z, t) 

= ]^da \dzA{s,z,t)\^ = 0 

since |9zA(s,z,t)p = by definition. The equations (4.15) and (4.16) imply 9sA(s,0,t) • 
5^A(5,0,t) = 0, which yields 

(4.17) 


^ 5 A(s, 2 ;, t) • dzA{s^ 2 ;, t) = 0 . 
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For simplification of the following calculations, we denote the variables s and z by si and S 2 
respectively. Given the arguments above, the metric tensor with respect to A is given by 

^11 — 9ss — dgA • 

gi2 = 921 = 9sz = 9zs = • dzA = 0 and 

922 = 9zz = dzA ■ dzA = 


The matrix 


f 911 

V^2i 922 J 


is thus diagonal and as usual we denote the entries of its inverse G~^ by 

For a scalar function /i(x, t) = h{s{t^ x)^z{t^ x)^t) on N ^ the surface gradient on C F is hence 

expressed by 


(4.18) Vvh = Y. g^^ds.hdsjA = g^^d,,hd,,A + -^ds,hds,A = + -d,hd,A 

i,j=i ^ ^ 


where is the surface gradient on = {A{s^ z^t)\s G U} for a fixed z G [0,1]. Similarly, 

(4.19) Vr • d — • a H —dzd • dzA 

for some vector valued function a{x^t) = a(s(x, t), 2 ;(x, t), t). 

In analogy to the appendix in [2], we calculate the Laplace-Beltrami operator Ap in the new 
coordinates. Due to the properties of the parametrization A{s^z^t) which already lead to the 
diagonal structure of the metric tensor above, we find • dzA = 0 and V^^/i • dzzA = 0. Hence 


( 9 ^ 7 ,^) • d,A = • d,A + • d,,A = d, {V-yJi ■ = 0 


and substituting (4.18) in (4.19) we thus compute 


Ap/i —A^^/i 


) (V,.9.S) . 


d,A + ■ d,A 


(4.20) 


+ I (V9.7e^) • d,A + ^d,Ji • 4A + ^dJid,,A ■ d,A 

=A^^/i H— dzh\/ • dzA H—• dzA 


where we have used the identities above. Since 

1 


dssA • dsA = -ds {dsA • dsA) = 0 , 


the curvature vector dssA of ^{t) (seen as a curve in M^) is an element in span (5^A, z/p) where z/p 
denotes the direction normal to the surface F. The geodesic curvature hig of ^{t) in F is therefore 
given by 

Kg = dssA • dzA = ds {dsA • dzA) - dsA • dszA = • a,A. 

As in [2], one can derive 

le ' h — ^ ^ ’ h 

A^^/i = A^/i + TZs 
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where TZ^ is of higher order in Thus (4.20) reads 


(4.21) 


Arh = Aj^h - ^Kgdzh + ^dzzh ■ dzA 

= Ajh - ^Kgdzh + ^dzzh ■ dzA + TZg. 


4.3. Asymptotic analysis: Inner expansion. We assume now that the functions in 0- 

(0 have inner expansions on N with respect to the new variables of the form 

oo 

fs{x,t) = F{z,s,t;s) = '^s''Fk{z,s,t) 

k=o 

where again ... etc. Accordingly, the inner expansions for 0, /i, 9) will be denoted 

by (y,$,M,©). 

In order for the inner and outer expansions to be consistent which each other, we prescribe the 
following matching conditions as 2 ; ^ ±oc (again we use F as an placeholder for (F, $, M, 0)) 

(4.22) Fo(t,s,±oo) ~ fo{x,t), 

(4.23) dzFo{t^ s, ±oc) ^ 0, 

(4.24) dzFi{t,s,±oo) ^Vrfo{x,t) 

where {x,t) = A(0, s,t) and f^{x,t) = lims^o fo{exp^{F5u-y),t). We refer the reader to [HI El] 
and m for a derivation of these matching conditions. 

We plug these asymptotic expansions in the equations and use the results from 

Remark |4.2| where necessary. Again we collect all terms with the same order in 6:. As we assume 
that the limes 6: ^ 0 exists, we require that the terms associated with the leading orders in 6: 
cancel out. The terms of order in equation ( |1.2[ ) hence yield the differential equation 

-5..$o + f^'(^o )=0 

for $0- Together with the matching conditions and the condition $o(0) = 0 we deduce 


(4.25) 


$o(^) = tanh( 2 ;) and 


-\dz^oWz) = W{^o){z). 


Looking at the conditions imposed by terms of order equation g yields OzzMq = 0 and 
integrating this equation from —00 to 2 ; implies that Mq is independent from 2 ; if we take the 
matching condition (4.23) into account. Thus 

Mq{z = +oc) = Mo{z — — 00 ), 


which in turn gives (4.9). 

In a similar way, we obtain equation (4.10). The terms of order in (1.7) indeed imply 
dzzQo — 0 and integrating in 2 : yields 


(4.26) 


dzQo = 0 


by the matching conditions. Again we deduce ©o(^ = + 00 ) = ©o(^ = — 00 ) and (4.22) yields 

Let us observe for later use that [0]^ = 0 directly implies 
(4.27) Hi = 1 

as we already saw that [(/:?] 1 = 2. 

A second observation is motivated by the study of all terms of order 0(1) in (1.3). We can 
deduce 

00 = I (2^0 - 1 - $ 0 ) 
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and since = 0 by (4.26) we therefore obtain 

(4.28) 2d;,Vo{z,s,t) = d:,^o{z). 

Substituting the inner expansions in ( |1.6[ ) also yields terms of order The resulting equation 
reads 

= dM. 


Equation (4.11) is then obtained from the matching conditions by integrating in 2 ;. 

Next, we study the terms of order 0(1) in (1.2). Similar to the studies on the related Cahn- 
Hilliard equation we obtain 

Mo = -^{2Vo-l- $ 0 ) 

= + Kgd.^o + W""($o)$i - ^00 

and apply the following solvability condition for $1 derived in [H Lemma 2.2]. 

Lemma 4.3 (|1]). Let A(z) be a hounded funetion on —00 < 2 ; < oc. Then the problem 


d,,^ + W"{Fo{z))^ = A{z) 


m = 0 


2 : G 
G 


has a solution if and only if 


A{z)Fq{z) dz = 0. 


It is indeed easy to see that this condition is necessary if one multiplies the equation by Fq and 
uses integration by parts. The assertion that the condition is also sufficient can be derived from 
the method of variation of constants, details are given in [3]. 

Lemma |4.3| directly yields 

1 

2Mo = Cong - - Qodz^o dz 

^ J —00 


where cq is given by 

/ oo 

dz. 

-OO 

Given the fact that ©0 is independent of 2 ; (see ( |4.26| )) we deduce 

/ oo roo 

Ood.^Q dz = eo / dz = 200 . 

OO J —00 

We thus conclude 

2 /i + 0 = 

as we claimed in (4.8). 

In order to conclude our analysis, we collect all terms of order from equation (1.7). We 
thus have 

and another integration with respect to 2 ; allows us to deduce 

-V [u]l = [Vr0]l • u.y 

from the matching conditions. Our observation (4.27|) on [n]l above hence implies (4.12). 
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4.4. Free energy inequalities and conservation properties. The discussion of the thermo¬ 
dynamical background in Section can be extended to the sharp interface model (4.2)-(4.12) 
derived above if one considers the surface free energy 


(4.29) 


'2 dn^ 


J='{(p,v,'y) ■=C{) J 1 ^ - 1 - U" 

, / 1 ^ + /+ 


= Co 


and the bulk free energy 
(4.30) 


J^5(r) '= - dx. 

2 J B 


In particular, choosing these energies the energy inequality derived in Lemma [2.1| for the diffuse 
interface model from the second law of thermodynamics has a counterpart in the sharp interface 
model, i.e. 

^ + < j^q{6- u) dn^. 

As the following calculations show, this is mostly due to the relations specified in the equations 


= / BgV dn\ 


(4.12) on 7 . They come into play since 

We refer the reader to Section 5.4 in m for a derivation of this identity. 

Let us now calculate the time derivative of (4.29) in order to derive the desired energy inequal¬ 
ity. Reynolds’ transport theorem and (4.31) imply 

d 
dt 


(4.32) 


= co J KgV diX + ~ (^J dn"^ ^ 


y2 ^ 


= co J KgV dn^ + ^J BdtO dn^ + ^J BdtO dH^ 


We have used equation ( 4.10[ ) to derive the last equality. Because of (4.8), the first integral 
Co KgV dlX coincides with y ( 2 /x + B)V dT-0-. Using (4.11) and (4.12) yields 


(4.33) 


Co [ KgV dn^ = - f n [Vr/r]l ■ dlX - f B [Vr0]l • dH^ 

J -J J ^ J ^ 

= — /iVrA • ^7 dT-L^ + / /rVr/i~ • dBX. 

J ^ J ^ 

- f B\7rB+ ■ dn^ + f 0Vr0“ • dU^. 

J 'y J y 


'7 *^7 

By the divergence theorem for tangential vector fields on manifolds, we see that all integrals 
equate to integrals over T^ and T” respectively. That is, we find 


/ /iVr/i^ • dH^ = — / div (/iVr/i) dT-L^ 

Jy Jr+ 


/ |Vr/ip dn^ - / /iAr/i dn^ 
/f+ Jt+ 


L 


|Vr/rp dn^ 
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The last equation holds because of (4.5). Keeping in mind that the orientation of 7 seen as the 
boundary of differs from the orientation of 7 seen as the boundary of T”, similar calculations 
lead to 


J MVr/x" ■u-ydy} = -J iVrMl^ 

For the other integrals in (4.33), we use ( |4.6| ) and (|^ to calculate 

- [ 6»Vr0+ ■v-ydy} = - [ |Vr6»|^ 

J'y JT+ 


dn^ 


OArO dn^ 


= - / iVr^r dn^ - e {dtv - q) dW 

Jt+ Jt+ 

= - [ |Vr0p dn^-^ [ OdtB dU^ + f Bq dU 

and the analogue result 

J BVtB- -p^dn^ = - j |Vr0| 

Plugging these findings in ( |4.33D yields 


dn^ [ BdtB dn^ + [ Bq dn^. 
4 Jr- Jv- 


Co 


J KgV dn^ = - J I Vr/i|^ dn^ - J |Vr0p dn^-^J BdtB dU^ 


+ J dq dn^ 


and by (4.32) we deduce 

= - J dn"^ - J iVrBl"^ dn^ + J Bq dU^ 

// 


<Jdq dH\ 

Similarly to the derivation of Lemma |2.1[ we also find 




dn^. 


This leads to the free energy inequality for the sharp interface model stated in |4.1 
(4.34) 


d 


— [T{(p,v,'y) + Tb{u)]< J^q{B-u)dn. 


Furthermore, equation (4.2) and Reynolds’ transport theorem allow us to deduce 


(4.35) 


4 /9<i«^=4 


dt 


dt 


T+(7 


IdU^-^ 
dt 


1 dn^ 


'r-(t) 


f v-u^dn^- [ v-u^ dn^^ 0 . 


We remark that the signs in the above equation depend on the values of (p in F+(t) and F“(t) 
respectively as well as the fact that 7 (t) is seen as the boundary of F+(t) for the first integral 
and as the boundary of F“(t) for the second integral. 


For energy densities chosen according to ( |4.29| ) and JT30|), equations ( |4.3[ ) and ( |4.4[ ) correspond 
to the mass balance equation for cytosolic cholesterol (| 2 . 11 ), while the mass balance equation ( 2 . 2 ) 
for the membrane bound cholesterol corresponds to (4.6). We use these equations to deduce 


(4.36) 


d 

dt 


/ u dx -\- V dH^ 
Ub Jr 


= 0 . 


Equations (4.35) and (4.36) show that the conservation laws derived in Lemma 2.1 hold for the 
sharp interface model as well. Equation (4.34) is the analogue to the general energy inequality in 
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Lemma |2.1| for the sharp interface model, if one chooses the constitutive relations in such a way 


that the resulting free energy densities lead to (4.29) and (4.30). 


5. Numerical simulations 


In this section we present numerical results for the reduced non-local system (3.8)-(3.11) and 


for the fully coupled bulk-surface model (1.2)-(1.8). For the first, we propose a discretization 


which is semi-implicit in time and which is based on surface finite elements in space OdS]. The 
details of the discretization are shown in Section IhTl We validate the numerical method in Section 
|5.2| with two benchmark tests and subsequently, in Section |5.3| we present simulation results 


showing properties of the model. Moreover, in Section [5M| we present a numerical simulation for 
a diffuse interface approximation of the fully coupled bulk-surface model For a large 

bulk diffusion coefficient, the results show a good agreement with those obtained for the reduced 


model in Section [573| This further justifies to use numerical simulations for the reduced model in 
order to study the properties of the original model. 

Cahn-Hilliard systems on stationary surfaces have been numerically investigated by parametric 
finite elements [l3] , level set techniques [23] and with diffuse interface methods [l3] mostly applied 
for the simulation of spinodal decomposition and coarsening scenarios. For corresponding results 
on moving surfaces we refer to [161 ES] (sharp interface approach) and [52] (diffuse interface 
approach), respectively. Recently, a Cahn-Hilliard like system has been numerically studied for 
the simulation of the influence of membrane proteins on phase separation and coarsening on cell 
membranes 


5.1. Surface FEM Discretization. For the discretization of the reduced system (3.8)-(3.11) 


with q given by (1.8) we develop a scheme similar to the one in m, where a related non-local 
reaction-diffusion system has been numerically treated. To discretize in time, we introduce steps 
Tm^ m = l,...,Mr. Denoting by the time discrete solution at time 

m 

tm — R, the time discrete weak formulation then reads 
1=1 


(5.1) 

(5.2) 

(5.3) 


f + f • VrV- = f —V-, 

Jr 'Un+i Jr Jr ^m+i 

- j + £ J • Vrr? + t J ? J 

+ i i y (W^"(<^("»))^("*) - -]fv, 

[ + 1 f • VrC - I / • VrC 

Jr Un+i ^ JF ^ Jr 

+ f = f —C + f 

Jr Jr ^m+i Jr 


for all G iL^(F), where the non-local term 


(5.4) 


(m) _ 

\B\ 


u' ' = 


M 


L 


(m) 


is treated fully explicitly. 

For the spatial discretization, we introduce a triangulation F/^ of the surface F, and we use linear 


surface finite elements to obtain from (5.1)-(5.4) a linear system of equations for the coefficients 
of the discrete solutions ^(m+i)^ ^6^+1) respect to the standard Lagrange basis. The 

resulting linear system is solved by a direct solver or a stabilized bi-conjugate gradient (BiCGStab) 
method depending on the number of unknowns. Furthermore, we apply a simple time adaptation 
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strategy, where time steps r^+i are chosen (within bounds) inversely proportional to 


max 

xeVh 


see e.g. 


. The above numerical scheme is implemented in the adaptive FEM toolbox AMDiS 


5.2. Benchmark tests. In this section we consider a triangulation F/^ = of the unit sphere 
5^, and we use the following set of parameters: 


(5.5) 


Cl = C 2 = 500, e = 0 . 02 , M = 5\B\ = S = 0 . 02 . 

O 


5.2.1. Analytic solution for a prescribed right hand side. In order to verify the validity of the 
above numerical scheme, we define right hand sides for (3.8) and (3.10) such that an analytic 
solution for the resulting system can be given. To be more precise, we let 

' 'd{x) + (3 — t^ 


(5.6) 


^p{x, t) := tanh 


V2e 


for X G 3“^ and t G [0,T). Moreover, (3 is an angle related to the initial value (p{x,0), and by 
i} = i3{x) we denote the angle given by 

'd{x) := arccos(a; 3 ) G [0 , tt] 

for x = {x\,X 2 .,xz) G S‘^. Furthermore, we let v{x,t) ■= |(1 + such that 6 = 0 holds. 

Then we define Fi, F 2 : 3“^ x [0,T] ^ M such that 

(5.7) dfifi — Ar/i = Fi 

(5.8) fj, =—sArif + e~^W'{(p) 

(5.9) dfV - q{u,v) = F 2 
hold with u = u{t) given by 


on X (0,r], 
on X (0,r], 
on 3^ X (0,r] 


(5.10) 


u= — \M 
47r 


-/ 


In this way we obtain expressions for Fi and F 2 , which are incorporated in the numerical scheme 
described in Section 5.1 Finally, we use the adopted scheme to find approximate numerical 
solutions iphi fo (5.7)-(5.10) for initial conditions for Vh given by corresponding values 

determined by (5.6) and v — 4(1 + (^). Note that the expression for F 2 involves an integral of (p 
which is numerically approximated in polar coordinates in every time step. 

With the usual Lagrange interpolation operator we define the relative errors 

,, \\hFi-,t) - Fhi-,t)\\L^{sF 

eooit) :=- 


ei{t) := 




of (ph{-,t) in the L°°(5'|)-, ^(^D-norms, respectively, and we consider their values for simulations 

with T = P = —j and different spatial resolutions in Table The results indicate a linear 
dependence of the L°°{3‘fP- and i7^(5'|)-errors, respectively, on the grid-size, as expected for 
elliptic problems. As a further illustration of this benchmark, in Figure [T| we see contour plots 
of <y 9 (-, 0 ) and (p{-,T), respectively. We remark that in the following we use a grid with 98306 
vertices for all simulations with spherical geometry. 
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number of vertices 

6146 

24578 

98306 

393218 

eoo(U 

0.245005 

0.0284765 

0.0107314 

0.00568034 

ei{T) 

0.2284131488 

0.0747088125 

0.0407739003 

0.0211030453 


Table 1. 


Relative errors for simulations of (5.7)-(5.10) with different resolutions. 



= 0) 



phi 


.OOOe+00 

;o.8 


I 

[ 0.4 
tO 
--0.4 


I 


■ 0.8 

■l.OOOe+00 


Figure 1. Benchmark simulation with prescribed right hand side in (5.7)-(5.9). Contour plots 


of initial values of analytic solution (5.6) and of values at time T. 


5.2.2. Validation of f v. We return to (3.8)-(3.11) and remark that one easily verifies from (3.12) 
that 




fulfills the ODE 
(5.11) 


Miri / iri M\ Cl 2 

. = C1^-(^C1^ + C2 + Cl^j.+ ^. . 


We compare solutions of the above ODE with 




obtained by numerical solution of ( |5.1[ )-(5.4), where we have chosen the initial conditions 

(^(., 0) = -0.25 + 7^, 0) = 0.25. 

Thereby, TZ :Th ^ [—0.001, 0.001] provides a perturbation given by an (irregular and nonperiodic) 
oscillation around zero. In Figure^ = f ^2 Vh obtained by numerical simulation of ( |5.![ )-([ 




is compared with a numerical solution of the ODE (5.11) with initial condition 


z{0) = f p(.,0) 

Js‘^ 

obtained with Maple^^[T]. The excellent agreement further justifies the chosen numerical scheme. 

5.3. Simulation results. 


5.3.1. Variation of parameters. We use the basic set of parameters given in (5.5). Furthermore, we 


use the initial conditions (/^(•, 0) = (f+TZ^ t(*, 0) = |, where (f = —0.5 and 7Z :Th ^ [—0.001, 0.001] 
again denotes a perturbation by an (irregular and nonperiodic) oscillation around zero. In Figure 
1^ we see corresponding results, where contour plots of the discrete solutions (first row) 

and Vhi’^t) (second row) are displayed for various times t. The evolution shows a spinodal 
decomposition with subsequent interrupted coarsening scenario. 

The geometric shape of the particles can drastically change, if one changes the values of (/;(•, 0). 
In a second example, we have used (^(-,0) =0 + 7^, while all remaining parameters have not been 
changed. The corresponding numerical results can be seen in Figure]^ 
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Figure 2. Comparison of a numerical solution of the ODE (5.11) with f ^2 Vh obtained by 
simulation of (|3.8|-(|3.11|). 






0.1177) 402.4228) 


| 1.024e-(00 
0.8 

0.4 


-0 


= 0) 


t ^ 0.0036) 


^ 0.0107) 




- 1.485e-(00 

I 

■^0.75 


I 


0.5 

2.500e-01 


Vh{‘,t = 0) 


Figure 3. Numerical results for (p(-,0) = —0.5 + IZ. First row: contour plots of for 

several choices of times t; second row: corresponding contour plots of u(-,t). 


In the following, we compare almost stationary states for simulations with varying parameters, 
in order to illustrate the properties of the model and its solutions. We use the basic set of 
parameters for the results shown in Figure and modify one parameter. Let 


(p 



Varying (p. We change (p by changing the initial condition for or rather p. In Figure one 
observes the almost stationary states from the previous two examples with p = 0 and p = —0.5 
and examples with intermediate values, where one can see a crossover from circular lipid rafts to 
stripe like patterns. For p = —0.75 the system allows for a constant stationary state with order 
parameter away from the classical phases p G { — 1,1}. 

Varying 6. For p = —0.5 we consider almost stationary discrete order parameter p^ for different 
values for the parameter 6. In Figure for large 6 the almost stationary p^ has one lipid raft, 
as one would expect for classical Cahn-Hilliard system on the sphere. For decreasing but still 
positive values of 6 we expect to approach stationary solution of Ohta-Kawasaki based dynamics, 
as shown in Section |3.4[ For a more detailed comparison with almost stationary states of Ohta- 


Kawasaki based dynamics we refer to the subsequent Section 5.3.2 
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= (fh{‘,t^ 0.003S) (f h{-,t ^ 0.0125) ^ 0.1035) (f h{-,t ^ 379.54:65) 



Vh{-,t = 0) Vh{-,t ^ 0.003S) Vh{‘,t ^ 0.0125) 0.1035) 379.5465) 


Figure 4. Numerical results for (y^(-, 0) = 0 + 77.. First row: contour plots of 9?(-, t) for several 
choices of times t; second row: corresponding contour plots of v{-,t). 




6 = 0.3 


M 

> , 


U' 

.f 


L 


J 


6 = 0.1 



6 = 0.02 



6 = 0.002 


I l.OlSe+OO 
0.8 

0.4 


j~-0.4 

E-0.8 

■^-9.731e-01 


Figure 6 . Almost stationary discrete solutions (fh for different values of 6. 


Varying ci. Returning to the previous standard parameters (p = —0.5 and 6 = 0.02 we investigate 
the influence of different values of ci on the almost stationary discrete states. Increasing ci 
corresponds to increasing the number of lipid rafts as shown in Figure 

Varying C2. With our standard choice ci = 500 we obtain for variation of values for C2 similar 
results as for the previous examples (varying ci). In Figure]^ one observes increasing number of 
rafts of decreasing sizes for increasing C2. 

From the analysis in Section 3^ we conclude that increasing ci or C2 corresponds to increasing 
the non-local energy contribution in the Ohta-Kawasaki functional (3.22), see also (3.23), which 
is expected to describe the dynamics in the limit 5^0. This fact provides an explanation of the 
increasing number of particles for increasing ci, C2. 
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Cl 


5 

Figure 7. 


Cl = 100 


Cl = 500 


Cl = 2000 


Almost stationary discrete solutions (fh for different values of ci. 


1 1.0090+00 


-|0 


^-0.4 

*-9.6390-01 



C2 = 5 


C2 = 100 


C2 = 500 


C2 = 2000 


1 1.0110+00 


■fO 


^-0.4 

*-9.5130-01 


Figure 8. Almost stationary discrete solutions (fh for different values of C 2 . 


5.3.2. Comparison with Ohta-Kawasaki model. For the result in Figure]^ we have investigated 
the following scenario. We run a simulation of (3.8)-( [3dT ) with parameter values as in (5.5) 
except 6 = 0.0001 towards an almost stationary state (see Figure]^ left). With the resulting 
discrete order parameter we continue a simulation with Ohta-Kawasaki-based dynamics with 
parameters according to (3.23) and (3.24) again towards an almost stationary state (see Figure 
middle). The difference between these two stationary results for is displayed in Figure 


(right). These results confirm the analytic results of Section 3.4 


A 

O 


^ c 

) < 

c 




w 



O 



) < 

C 


L. 





ll -0.0004 
*-7.9960-04 


Figure 9. Almost stationary (fh obtained from a simulation of the reduced system (left) 
and from subsequent Ohta-Kawasaki-based dynamics (middle), difference between the previous 
numerical solutions (right). 


5.3.3. Non-spherical membrane shape. In a further example, we show in Figure 10 results with 


the standard parameter set (5.5) for the case that F is not a sphere. One obtains similar patterns 
as for the sphere. However, this configuration appears to be more stable than the previous 
spherical ones, where sometimes very slow arrangements could take place. This behavior has not 
been observed for this geometry. We remark that for this simulation, we have chosen a similar 
resolution as for the previous examples with spherical geometry. 

5.3.4. Comparison with energy deereasing dynamies. Here, we study a choice of q leading to an 
energy decreasing evolution. From the analysis in Section |3.3| we expect that solutions of the 


system (3.8)-(3.11) generically converge to configurations with one lipid phase concentrated in a 


single geodesic ball. For this purpose, we choose 

q — —c{0 — u) 














SURFACE CAHN-HILLIARD MODELING RAFT FORMATION 


27 



t = 0) t ^ 0.0053) t ^ 1.4133) t ^ 55.3218) 


Figure 10. Numerical results for a non-spherical surface and 7^(-,0) = —0.5 + IZ. Contour 
plots of for several choices of times t. 


with c = 500, corresponding to (2.17) and fb{u) = in (2.15). The results in this case are 
displayed in Figure [m and illustrate the evolution towards an almost stationary state with a 
single spherical raft particle. 



= 0) 



(ph(-, t ^ 0.0030) 



0.0103) (f h{-,t ^ 1.0166) ^ 2.0562) 


■ 1.038e+00 

Fo.8 

||§0.4 



Figure 11. Numerical results for a choice of q leading to an energy decreasing evolution and 
0) = —0.5 + 7Z. Contour plots of t) for several choices of times t. 


5.4. Simulation of the full model. For the simulation of the coupled bulk-surface model 


(1.2)-(1.8), we propose a diffuse interface approximation based on the diffuse approaches for the 


treatment of PDF’s on and inside closed surfaces in [43] and [33|, respectively. For a related phase- 
field description of a coupled bulk-surface reaction-diffusion system, we refer to m- Moreover, 
a further diffuse interface description of a coupled bulk-surface system has been given in m, 
and in [3] a diffuse interface approach of a linear coupled elliptic PDF system has been analyzed 
with respect well-posedness of the diffuse system and convergence towards its sharp interface 
counterpart. 

We choose a domain D C containing B and provide a phase-field approximation of 0 - 
0 by 


(5.12) 

(5.13) 

(5.14) 

(5.15) 

where the phase-field function ip : fl 


(5.16) 


ip{x) := 


in D X (0, T], 
in D X (0, T], 
in D X (0, T], 

in D X (0, T], 

describing the geometry, that is B and P, is given by 
3r{x) 


^dtu = DV • {'ipVu) — Eg ^h{^)q{u^ v) 

= V • (5(^)V/i) 

b{ip)ii = -eV • {b{ip)Vip) + e-^b{ip)W'{ip) - 5-^b{ip){2v - 1 - if) 

b{ip)dtv = |v • {b{iP)Vv) - |v • {b{iP)Vip) + b{ip)q{u,v) 

0 0 


1 — tanh 


for a (small) real number > 0 and a signed distance r to F being negative in B. Note that ?/) 
is obtained by smearing out the characteristic function xb of S on a length proportional to Eg. 
Moreover, b = b{^) := is small outside the diffuse interface region. 

We use a semi-implicit adaptive FEM discretization (see [12]) and choose D = (—2,2)^ with 
all solutions assumed D-periodic, Eg = |, D = 100. All otherwise degenerate mobility functions 


appearing in second order terms in (5.12)-(5.15) are regularized by addition of a parameter 
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6r = 10“^. The resulting linear system is solved by a stabilized bi-conjugate gradient (BiCGStab) 
solver. The numerical scheme is implemented in AMDiS j50j . Moreover, we apply the parameters 
in (5.5) and initial conditions as in Section 5.3 and additionally r(-, 0) = 4.25 corresponding to 

plots of foi* various times where 


M in 


(5.5). In Figure 


12 


we see 


is obtained by replacing r by adiscrete signed distance in (|5.16 ). The results are in good 
qualitative agreement with the results in Section [5^ (see Figure and hence further justify the 
reduced model formally obtained in the limit D 


oo. 



= 0) 
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0.0126) 
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Figure 12. Numerical results for diffuse interface approximation of the full model (1.2)-(1.8). 
Contour plots of i} for several choices of times t. 


6. CONGLUSIONS 

We have presented a model for lipid raft formation that extends in particular the model by 
Gomez, Sagues and Reigada [22]. The key new aspect in our work is an explicit account for 
the cholesterol dynamics in the cytosol, which leads to a complex system of partial differential 
equations both in the bulk and on the cell membrane. These are coupled by an outffow condition 
for the cytosolic cholesterol and a source term in the surface membrane equation, determined by 
a constitutive relation. These considerations lead to an interesting and complex mathematical 
model. The surface equations combine a Cahn-Hilliard type evolution equation for an order 
parameter and an equation for the surface cholesterol. The latter is a diffusion equation including 
a cross-diffusion term and an additional (nonlocal) term. In the bulk we have a diffusion equation 
with a Robin-type boundary condition that couples both systems. 

We have shown that our model can be derived from thermodynamical conservation laws and 
free energy inequalities for bulk and surface processes. Depending on the specific choice of the 
exchange term we obtain a total free energy decrease or not, which can be seen as a distinction 
between equilibrium (closed system) and non-equilibrium (open system) type models. The anal¬ 
ysis of both classes reveals striking differences in the behavior: whereas equilibrium-type models 
only support the formation of macrodomains and connected phases, a prototypical choice of a 
non-equilibrium model leads to the formation of raft-like structures. 

These findings are the result of both a (formal) qualitative mathematical analysis and careful 
numerical simulations. More precisely, in a certain parameter regime (large interaction strength 
between lipid and cholesterol) we have demonstrated that stationary states of the prototypical 
non-equilibrium model are close to stationary states of the Otha-Kawasaki model, whereas in 
the case of equilibrium models stationary states coincide with stationary states for a surface 
Cahn-Hilliard equation. 

For a better understanding of the (complex) model asymptotic reductions are instrumental. 
We in particular justify the sharp interface limit of the model, that is the reduction when the 
parameter s (related to the width of transition layers between the distinct lipid phases) tends to 
zero. Moreover, we have derived a simplified model by taking the cytosolic diffusion (typically 
much larger than the lateral membrane diffusion) to infinity. This reduction leads to a non-local 
model that only includes variables on the surface membrane. 
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Numerical simulation show that, depending on the parameter choices, uniformly distributed 
circular microdomains or stripe pattern that stretch out over the membrane emerge. A justifica¬ 
tion of the numerical scheme and implementation is presented and the dependence on the various 
parameters is analyzed. Further evidence for the connection with the Ohta-Kawasaki dynamics 
(with specific parameters that are in functional relation to the parameters of our model) is given. 

Our results support the belief that non-equilibrium processes are indeed essential for mi¬ 
crodomain formation. On the other hand, emerging structures are—in the long-time behavior— 
very regular and uniform, in contrast to the picture of a persisting redistribution of location and 
sizes of rafts observed in experiments. Thermal fluctuations and/or the interaction with proteins 
on the membrane have to be taken into account to obtain this more complex behavior. Our 
contribution presents a solid basis for such extensions. 
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